output_version <- '2024-07-21' #
## fishnet
fishnet_lake_chad_ntl_adm0_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01')
fishnet_adm0_only <- read_sf(dirname(fishnet_lake_chad_ntl_adm0_x_d01_shp),basename(fishnet_lake_chad_ntl_adm0_x_d01_shp))
st_crs(fishnet_adm0_only) <- "EPSG:4326"


##
## Source: Gridded Livestock of the World (GLW 3) database, 
## reflecting the most recently compiled and 
## harmonized subnational livestock distribution data for 2010
## For each cell (same cell, 50km, 100km), we need the values 
## of the different livestock types and ideally if there is 
## an aggregation in terms of TLU = total livestock units = 
## using the coefficients for Africa from annex 1 in this document: 
## https://www.fao.org/4/i2294e/i2294e00.pdf
##

# livestock - cattle
cattle_fn <- file.path(wkdir,sprintf('local/gis_data/001_farming/GLW/livestock_%s','6_Ct_2010_Aw.tif'))
xfun::dir_create(dirname(cattle_fn))

if(!file.exists(cattle_fn)){
  url_cattle <- 'https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/5U8MWI/T7QSGJ'
  response <- httr::GET(url_cattle)
  if (httr::status_code(response) == 200) {
    # Save the file
    writeBin(content(response, "raw"), cattle_fn)
    cat("File downloaded successfully!")
  } else {
    cat("Error downloading file:", http_status(response)$message)
  }
  lstk_ct <- terra::rast(cattle_fn)
} else {
  lstk_ct <- terra::rast(cattle_fn)
}

# livestock - chicken
chicken_fn <- file.path(wkdir,sprintf('local/gis_data/001_farming/GLW/livestock_%s','6_Ch_2010_Aw.tif'))

if(!file.exists(chicken_fn)){
  url_chicken <- 'https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/SUFASB/AP1LHN'
  response <- httr::GET(url_chicken)
  if (httr::status_code(response) == 200) {
    # Save the file
    writeBin(content(response, "raw"), chicken_fn)
    cat("File downloaded successfully!")
  } else {
    cat("Error downloading file:", http_status(response)$message)
  }
  lstk_ch <- terra::rast(chicken_fn)
} else {
  lstk_ch <- terra::rast(chicken_fn)
}

 	

# livestock - horses
horses_fn <- file.path(wkdir,sprintf('local/gis_data/001_farming/GLW/livestock_%s','6_Ho_2010_Aw.tif'))

if(!file.exists(horses_fn)){
  url_horses <- 'https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/7Q52MV/UAWH3Z'
  response <- httr::GET(url_horses)
  if (httr::status_code(response) == 200) {
    # Save the file
    writeBin(content(response, "raw"), horses_fn)
    cat("File downloaded successfully!")
  } else {
    cat("Error downloading file:", http_status(response)$message)
  }
  lstk_ho <- terra::rast(horses_fn)
} else {
  lstk_ho <- terra::rast(horses_fn)
}


# livestock - goats
goats_fn <- file.path(wkdir,sprintf('local/gis_data/001_farming/GLW/livestock_%s','6_Gt_2010_Aw.tif'))

if(!file.exists(goats_fn)){
  url_goats <- 'https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/OCPH42/ZMVOOW'
  response <- httr::GET(url_goats)
  if (httr::status_code(response) == 200) {
    # Save the file
    writeBin(content(response, "raw"), goats_fn)
    cat("File downloaded successfully!")
  } else {
    cat("Error downloading file:", http_status(response)$message)
  }
  lstk_gt <- terra::rast(goats_fn)
} else {
  lstk_gt <- terra::rast(goats_fn)
}



# livestock - pig
pig_fn <- file.path(wkdir,sprintf('local/gis_data/001_farming/GLW/livestock_%s','6_Pg_2010_Aw.tif'))

if(!file.exists(pig_fn)){
  url_pigs <- 'https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/33N0JG/9LGGBS'
  response <- httr::GET(url_pigs)
  if (httr::status_code(response) == 200) {
    # Save the file
    writeBin(content(response, "raw"), pig_fn)
    cat("File downloaded successfully!")
  } else {
    cat("Error downloading file:", http_status(response)$message)
  }
  lstk_pg <- terra::rast(pig_fn)
} else {
  lstk_pg <- terra::rast(pig_fn)
}


# livestock - sheep
sheep_fn <- file.path(wkdir,sprintf('local/gis_data/001_farming/GLW/livestock_%s','6_Sh_2010_Aw.tif'))

if(!file.exists(sheep_fn)){
  url_sheep <- 'https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/BLWPZN/XDIRM4'
  response <- httr::GET(url_sheep)
  if (httr::status_code(response) == 200) {
    # Save the file
    writeBin(content(response, "raw"), sheep_fn)
    cat("File downloaded successfully!")
  } else {
    cat("Error downloading file:", http_status(response)$message)
  }
  lstk_sh <- terra::rast(sheep_fn)
} else {
  lstk_sh <- terra::rast(sheep_fn)
}


#For each cell (same cell, 50km, 100km), we need the area shares of the various categories below:

require(exactextractr)

sum_cover <- function(x){
  list(x %>%
         group_by(value) %>%
         dplyr::summarize(total_area = sum(coverage_area)) %>%
         dplyr::mutate(sh_landsuit = total_area/sum(total_area)))
}


##
## RADIUS BUFFER 0 (grid cell only) 50KM (FROM POINT), 100KM (FROM POINT)
##

livestock_files <- list.files(file.path(wkdir,'local/gis_data/001_farming'),
                              pattern='livestock', full.names = T) 

livestock_TLU_coefficients <- openxlsx::read.xlsx(file.path(wkdir,
                                                            'local/gis_data/001_farming',
                                                            'TLU_units.xlsx'))
## memory issues
n_objectid_list=unique(fishnet_adm0_only$OBJECTID)
total_chunks <- split(n_objectid_list, ceiling(seq_along(n_objectid_list)/1000))
length(total_chunks)

# process over livestock #
for(livestock_fn_i in livestock_files){
  print(livestock_fn_i)

  ## START PROCESSING ##
  tif_fn <- substr(basename(livestock_fn_i),11,nchar(basename(livestock_fn_i)))
  
  # livestock type from lookup table
  livestock_type <- livestock_TLU_coefficients[which(livestock_TLU_coefficients$livestock_aw_tif==tif_fn),"livestock_type"]
  
  # livestock TLU coefficient from lookup table
  livestock_TLU_value <- livestock_TLU_coefficients[which(livestock_TLU_coefficients$livestock_aw_tif==tif_fn),"tlu_parameter"]
  if(!is.na(livestock_TLU_value)){
    livestock_TLU_i <- livestock_TLU_value
  } else {
    # else default 0
    livestock_TLU_i <- 1  
  }
  
  ## memory issues
  library(doParallel)
  
  for(chunk_index in 1:length(total_chunks)) {
    print(chunk_index)
    livestock_chunk_out_fn <- file.path(wkdir,'proc_data',
                                       'livestock_chunks',sprintf('LIVESTOCK_%s_TLU_%s_of_%s_%s.csv',livestock_type,chunk_index,length(total_chunks),output_version))
    xfun::dir_create(dirname(livestock_chunk_out_fn))
    if(!file.exists(livestock_chunk_out_fn)){
      n_objectid_list_chunk <- as.numeric(unlist(total_chunks[chunk_index]))
      #n_objectid_list_chunk=1:2
      
      n_cores <- min(10,length(n_objectid_list_chunk))
  
      cl <- makeCluster(n_cores)
      registerDoParallel(cl)
      
      results <- foreach(objectid_j=n_objectid_list_chunk,
                         .export=c(ls()),
                         .packages = c("dplyr", "sf","exactextractr","spatialEco","terra")) %dopar% {
                           # Modulus operation
                           if(objectid_j %% 100==0) {
                             # Print on the screen some message
                             cat(paste0("iteration: ", objectid_j, "\n"))
                           }
                           
                           fishnet_lake_chad_ntl_adm0_x_d01_shp_fn <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01.shp')
                           fishnet_adm0_only_j = read_sf(fishnet_lake_chad_ntl_adm0_x_d01_shp_fn, query = sprintf('SELECT * from fishnet_lake_chad_ntl_adm0_x_d01 WHERE OBJECTID = %s',objectid_j))
                           st_crs(fishnet_adm0_only_j) <- "EPSG:4326"
                           
                           if(dim(fishnet_adm0_only_j)[1]>0){
                             
                             fishnet_adm0_only_buf100km_bbox <- round(st_bbox(spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=150000)),0)
                             ## crop raster to 100km buffer
                             livestock_raw_i <- terra::rast(livestock_fn_i)
                             
                             # apply scalar TLU Total Livestock Unit coefficient
                             livestock_tlu_i <- livestock_raw_i * livestock_TLU_i
                             
                             lstk_crop_j <- terra::crop(livestock_tlu_i,fishnet_adm0_only_buf100km_bbox,snap="out")
                             all_na <- all(is.na(values(lstk_crop_j)))
                             if(all_na==FALSE){
                               ## 0km
                               x <- exact_extract(lstk_crop_j, fishnet_adm0_only_j, 'sum')
                               
                               #merge the list elements into a df
                               livestock_TLU_0 <- as.data.frame(cbind(x, objectid = objectid_j)) %>%
                                 mutate(buffer_km = 0)
                               
                               rm(x)
                               # projected 50km
                               fishnet_adm0_only_buf50km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=50000)
                               x <- exact_extract(lstk_crop_j, fishnet_adm0_only_buf50km, fun = 'sum')
                               
                               #merge the list elements into a df
                               livestock_TLU_50 <- as.data.frame(cbind(x, objectid = objectid_j)) %>%
                                 mutate(buffer_km = 50)
                               
                               # projected 50km
                               fishnet_adm0_only_buf100km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=100000)
                               x <- exact_extract(lstk_crop_j, fishnet_adm0_only_buf100km, fun = "sum")
                               #merge the list elements into a df
                               livestock_TLU_100 <- as.data.frame(cbind(x, objectid = objectid_j)) %>%
                                 mutate(buffer_km = 100)
                               
                               out <- rbind(livestock_TLU_0,livestock_TLU_50,livestock_TLU_100)
                               names(out)[1] <- "livestock_TLU"
                               
                               # add objectid
                               out$objectid <- objectid_j
                               out$livestock_source <- tools::file_path_sans_ext(tif_fn)
                                 
                               # clean-up
                               rm(x,fishnet_adm0_only_buf100km,fishnet_adm0_only_buf50km)
                               
                               return(out)  
                             }
                           }
                           
                         }
      
      parallel::stopCluster(cl)
      
      system.time(livestock_dt <- data.table::rbindlist(results))
      livestock_df <- livestock_dt %>% as.data.frame() %>%
          dplyr::select(objectid,livestock_TLU,livestock_source,buffer_km) %>%
          arrange(objectid,buffer_km,livestock_source)
      names(livestock_df)
      summary(livestock_df)
      
      data.table::fwrite(livestock_df,livestock_chunk_out_fn)
      
      # clean-up
      rm(livestock_df,livestock_dt)
    }
      
    
  }
  
  # merge all chuncks
  livestock_df_files <- sort(list.files(dirname(livestock_chunk_out_fn),
                                        pattern=sprintf("LIVESTOCK_%s",livestock_type),
                                        full.names=T))
  tables <- lapply(livestock_df_files, function(z) read.csv(z))
  livestock_df_all <- data.table::rbindlist(tables)
  
  #
  livestock_out <- livestock_df_all %>%
    dplyr::select(objectid,livestock_TLU,livestock_source,buffer_km) %>%
    arrange(objectid,buffer_km,livestock_source)
  
  # out file  
  livestock_out_fn <- file.path(wkdir,'proc_data',sprintf('LIVESTOCK_TLU_%s_%s.dta',livestock_type, output_version))
  
  # Adding the label value
  attr(livestock_out$objectid, "label") <- "objectid of grid"
  attr(livestock_out$livestock_TLU, "label") <- 'Livestock Total Livestock Units (SSA) : GLW3 AW'
  attr(livestock_out$livestock_source, "label") <- 'Livestock type : GLW3 '
  attr(livestock_out$buffer_km, "label") <- 'Buffer'
  attr(livestock_out, "label") <- sprintf("livestock_production.R last run on %s",Sys.Date())
  head(livestock_out)
  summary(livestock_out)
  
  ## Save the data
  haven::write_dta(livestock_out, livestock_out_fn)
  
  ## END PROCESSING LIVESTOCK_i ##
}
